{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 实验七 线性方程组的迭代解法\n",
    "\n",
    "雅可比（Jacobi）迭代法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def jacobi_iter(A, b, x0=None, eps=1e-5, max_steps=5000, verbose=False):\n",
    "    \"\"\"雅可比（Jacobi）迭代法求解线性方程组:\n",
    "        A @ x = b\n",
    "    \n",
    "    Args:\n",
    "        A:  np_array_like, 系数矩阵\n",
    "        b:  np_array_like, 右端常数\n",
    "        x0: np_array_like, 迭代初值\n",
    "            default x0=None means using random values.\n",
    "        eps: float, 精度要求\n",
    "        max_steps: int, 最大迭代次数\n",
    "        verbose: bool, 如果计算成功，打印出结果及迭代次数\n",
    "        \n",
    "    Returns:\n",
    "        x: 方程组的解\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: A 和 b 形状不符合要求\n",
    "        Expection:  达到最大迭代次数，仍不满足精度\n",
    "    \"\"\"\n",
    "    A = np.array(A)\n",
    "    b = np.array(b)\n",
    "    \n",
    "    n, m = A.shape\n",
    "    if n != m or n != len(b):\n",
    "        raise ValueError(f\"Not match: A({n, m}) and b({len(b)},)\")\n",
    "        \n",
    "    if not x0:\n",
    "        x0 = np.random.random(n)\n",
    "        \n",
    "    # A =  D - L - U\n",
    "    \n",
    "    D = np.diag(np.diag(A))\n",
    "    L = - np.tril(A, -1)\n",
    "    U = - np.triu(A, 1)\n",
    "    \n",
    "    inv_D = np.linalg.pinv(D)\n",
    "    \n",
    "    B = inv_D @ (L + U)\n",
    "    f = inv_D @ b\n",
    "    \n",
    "    x = x0\n",
    "    i = 0\n",
    "    for i in range(int(max_steps)):\n",
    "        x_prev = np.array(x)  # deep copy\n",
    "        \n",
    "        x = B @ x + f\n",
    "        \n",
    "        if np.all(np.abs(x - x_prev) <= eps):  # 达到精度要求\n",
    "            break\n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\")\n",
    "    \n",
    "    if verbose:\n",
    "        print(f\"jacobi_iter get result x = {x} after {i} iterations.\")\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jacobi_iter get result x = [0.99999906 0.99999911 0.99999933] after 7 iterations.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.99999906, 0.99999911, 0.99999933])"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobi_iter([[9,-1,-1], [-1,10,-1],[-1,-1,15]], [7,8,13], verbose=True)\n",
    "# jacobi_iter([[9,-1,-1], [-1,10,-1],[-1,-1,15]], [7,8,13], x0=[99999, 99999, 99999], verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "高斯-赛德尔迭代"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_seidel_iter(A, b, x0=None, eps=1e-5, max_steps=5000, verbose=False):\n",
    "    \"\"\"高斯-赛德尔迭代（Gauss–Seidel method）求解线性方程组:\n",
    "        A @ x = b\n",
    "    \n",
    "    Args:\n",
    "        A:  np_array_like, 系数矩阵\n",
    "        b:  np_array_like, 右端常数\n",
    "        x0: np_array_like, 迭代初值\n",
    "            default x0=None means using random values.\n",
    "        eps: float, 精度要求\n",
    "        max_steps: int, 最大迭代次数\n",
    "        verbose: bool, 如果计算成功，打印出结果及迭代次数\n",
    "        \n",
    "    Returns:\n",
    "        x: 方程组的解\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: A 和 b 形状不符合要求\n",
    "        Expection:  达到最大迭代次数，仍不满足精度\n",
    "    \"\"\"\n",
    "    A = np.array(A)\n",
    "    b = np.array(b)\n",
    "    \n",
    "    n, m = A.shape\n",
    "    if n != m or n != len(b):\n",
    "        raise ValueError(f\"Not match: A({n, m}) and b({len(b)},)\")\n",
    "        \n",
    "    if not x0:\n",
    "        x0 = np.random.random(n)\n",
    "\n",
    "    # A =  D - L - U\n",
    "    \n",
    "    D = np.diag(np.diag(A))\n",
    "    L = - np.tril(A, -1)\n",
    "    U = - np.triu(A, 1)\n",
    "    \n",
    "    inv_DsL = np.linalg.pinv(D - L)\n",
    "    \n",
    "    B = inv_DsL @ U\n",
    "    f = inv_DsL @ b\n",
    "    \n",
    "    x = x0\n",
    "    i = 0\n",
    "    for i in range(int(max_steps)):\n",
    "        x_prev = np.array(x)  # deep copy\n",
    "        \n",
    "        inv_DsL = np.linalg.pinv(D - L)\n",
    "        x = B @ x + f\n",
    "        \n",
    "        if np.all(np.abs(x - x_prev) <= eps):  # 达到精度要求\n",
    "            break\n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\")\n",
    "    \n",
    "    \n",
    "    if verbose:\n",
    "        print(f\"gauss_seidel_iter get result x = {x} after {i} iterations.\")\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gauss_seidel_iter get result x = [1. 1. 1.] after 10 iterations.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1.])"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gauss_seidel_iter([[9,-1,-1], [-1,10,-1],[-1,-1,15]], [7,8,13], eps=1e-12, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SOR 迭代法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sor(A, b, x0=None, w=1, eps=1e-5, max_steps=5000, verbose=False):\n",
    "    \"\"\"逐次超松驰法（Successive over-relaxation, SOR）求解线性方程组:\n",
    "        A @ x = b\n",
    "    \n",
    "    Args:\n",
    "        A:  np_array_like, 系数矩阵\n",
    "        b:  np_array_like, 右端常数\n",
    "        x0: np_array_like, 迭代初值\n",
    "            default x0=None means using random values.\n",
    "        w:  float, 松弛因子 w > 0\n",
    "            default w=1，即 Gauss–Seidel 迭代\n",
    "        eps: float, 精度要求\n",
    "        max_steps: int, 最大迭代次数\n",
    "        verbose: bool, 如果计算成功，打印出结果及迭代次数\n",
    "        \n",
    "    Returns:\n",
    "        x: 方程组的解\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: w 小于等于 0\n",
    "        ValueError: A 和 b 形状不符合要求\n",
    "        Expection:  达到最大迭代次数，仍不满足精度\n",
    "    \"\"\"\n",
    "    if w <= 0:\n",
    "        raise ValueError(f\"unexpected w = {w} < 0\")\n",
    "    \n",
    "    A = np.array(A)\n",
    "    b = np.array(b)\n",
    "    \n",
    "    n, m = A.shape\n",
    "    if n != m or n != len(b):\n",
    "        raise ValueError(f\"Not match: A({n, m}) and b({len(b)},)\")\n",
    "        \n",
    "    if not x0:\n",
    "        x0 = np.random.random(n)\n",
    "\n",
    "    # A =  D - L - U\n",
    "    \n",
    "    D = np.diag(np.diag(A))\n",
    "    L = - np.tril(A, -1)\n",
    "    U = - np.triu(A, 1)\n",
    "    \n",
    "    _inv_DsWL = np.linalg.pinv(D - w * L)\n",
    "    \n",
    "    B = _inv_DsWL @ ((1-w) * D + w * U)\n",
    "    f = w * _inv_DsWL @ b\n",
    "    \n",
    "    x = x0\n",
    "    i = 0\n",
    "\n",
    "    for i in range(int(max_steps)):\n",
    "        x_prev = np.array(x)  # deep copy\n",
    "        \n",
    "        x = B @ x + f\n",
    "        \n",
    "        if np.all(np.abs(x - x_prev) <= eps):  # 达到精度要求\n",
    "            break\n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\")\n",
    "    \n",
    "    \n",
    "    if verbose:\n",
    "        print(f\"sor (w={w}) get result x = {x} after {i} iterations.\")\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sor (w=1.1) get result x = [1. 1. 1.] after 14 iterations.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1.])"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sor([[9,-1,-1], [-1,10,-1],[-1,-1,15]], [7,8,13], w=1.1, eps=1e-12, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jacobi_iter get result x = [-1. -1. -1. -1.] after 68 iterations.\n",
      "gauss_seidel_iter get result x = [-1. -1. -1. -1.] after 36 iterations.\n",
      "sor (w=1.3) get result x = [-1. -1. -1. -1.] after 20 iterations.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-1., -1., -1., -1.])"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[-4, 1, 1, 1], \n",
    "     [1, -4, 1, 1],\n",
    "     [1, 1, -4, 1],\n",
    "     [1, 1, 1, -4]]\n",
    "b = [1, 1, 1, 1]\n",
    "\n",
    "x0 = [0, 0, 0, 0]\n",
    "\n",
    "jacobi_iter(A, b, x0=x0, eps=1e-9, verbose=True)\n",
    "gauss_seidel_iter(A, b, x0=x0, eps=1e-9, verbose=True)\n",
    "sor(A, b, x0=x0, w=1.3, eps=1e-9, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "重复的代码有点多，面向对象封装一下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SimpleIteration(object):\n",
    "    \"\"\"求解线性方程组的迭代法:\n",
    "        A @ x = b\n",
    "    \n",
    "    调用 SimpleIteration 子类实例得到解, e.g. JacobiIteration(A, b)()\n",
    "    \n",
    "    \"\"\"\n",
    "    def __init__(self, A, b):\n",
    "        \"\"\"\n",
    "        A:  np_array_like, 系数矩阵\n",
    "        b:  np_array_like, 右端常数\n",
    "        \"\"\"\n",
    "        self.A = np.array(A)\n",
    "        self.b = np.array(b)\n",
    "\n",
    "        n, m = self.A.shape\n",
    "        if n != m or n != len(self.b):\n",
    "            raise ValueError(f\"Not match: A({n, m}) and b({len(b)},)\")\n",
    "        \n",
    "    @staticmethod\n",
    "    def _dlu(A):\n",
    "        \"\"\"分裂 A:\n",
    "            A = D - L - U\n",
    "\n",
    "        Args:\n",
    "            A: np.array\n",
    "\n",
    "        Returns:\n",
    "            (D, L, U)\n",
    "        \"\"\"\n",
    "        D = np.diag(np.diag(A))\n",
    "        L = - np.tril(A, -1)\n",
    "        U = - np.triu(A, 1)\n",
    "\n",
    "        return D, L, U\n",
    "    \n",
    "    def _B_f(self):\n",
    "        \"\"\"计算迭代 x = B @ x + f 的 B 和 f\n",
    "        \"\"\"\n",
    "        raise NotImplementedError('_B_f')\n",
    "                \n",
    "    def __call__(self, x0=None, eps=1e-5, max_steps=5000, verbose=False):\n",
    "        \"\"\"线性方程组「简单迭代法」的迭代过程\n",
    "            x = B @ x + f\n",
    "            \n",
    "        其中 B, f 调用 self._B_f() 得到\n",
    "        \n",
    "\n",
    "        Args:\n",
    "            x0: np_array_like, 迭代初值\n",
    "                default x0=None means using random values.\n",
    "            eps: float, 精度要求\n",
    "            max_steps: int, 最大迭代次数\n",
    "            verbose: bool, 如果计算成功，打印出结果及迭代次数\n",
    "\n",
    "        Returns:\n",
    "            x: 方程组的解\n",
    "\n",
    "        Raises:\n",
    "            ValueError: x0 形状不符合问题\n",
    "            Expection:  达到最大迭代次数，仍不满足精度\n",
    "        \"\"\"\n",
    "        if not x0:\n",
    "            x0 = np.random.random(self.A.shape[0])\n",
    "        \n",
    "        x0_shape = np.shape(x0)\n",
    "        if x0_shape[0] != self.A.shape[0]:\n",
    "            raise ValueError(f\"Not match: A({self.A.shape}) and x0({x0_shape})\")\n",
    "        \n",
    "        B, f = self._B_f()\n",
    "        \n",
    "        x = x0\n",
    "        i = 0\n",
    "        for i in range(int(max_steps)):\n",
    "            x_prev = np.array(x)  # deep copy\n",
    "\n",
    "            x = B @ x + f\n",
    "\n",
    "            if np.all(np.abs(x - x_prev) <= eps):  # 达到精度要求\n",
    "                break\n",
    "        else:\n",
    "            raise Exception(f\"{self.method_name()} cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\")\n",
    "\n",
    "            \n",
    "        if verbose:\n",
    "            print(f\"{self.method_name()} get result x = {x} after {i} iterations.\")\n",
    "            \n",
    "    def method_name(self):\n",
    "        return self.__class__.__name__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class JacobiIteration(SimpleIteration):\n",
    "    \"\"\"雅可比（Jacobi）迭代法求解线性方程组:\n",
    "        A @ x = b\n",
    "    \"\"\"\n",
    "    def _B_f(self):\n",
    "        D, L, U = self._dlu(self.A)\n",
    "        \n",
    "        inv_D = np.linalg.pinv(D)\n",
    "\n",
    "        B = inv_D @ (L + U)\n",
    "        f = inv_D @ self.b\n",
    "        \n",
    "        return B, f\n",
    "    \n",
    "class GaussSeidel(SimpleIteration):\n",
    "    \"\"\"高斯-赛德尔迭代（Gauss–Seidel method）求解线性方程组:\n",
    "        A @ x = b\n",
    "    \"\"\"\n",
    "    def _B_f(self):\n",
    "        D, L, U = self._dlu(self.A)\n",
    "        \n",
    "\n",
    "        inv_DsL = np.linalg.pinv(D - L)\n",
    "\n",
    "        B = inv_DsL @ U\n",
    "        f = inv_DsL @ self.b\n",
    "        \n",
    "        return B, f\n",
    "    \n",
    "class SOR(SimpleIteration):\n",
    "    \"\"\"逐次超松驰法（Successive over-relaxation, SOR）求解线性方程组:\n",
    "        A @ x = b\n",
    "    \"\"\"\n",
    "    def __init__(self, A, b, w=1):\n",
    "        \"\"\"\n",
    "        w:  float, 松弛因子 w > 0\n",
    "            default w=1，即 Gauss–Seidel 迭代\n",
    "        \"\"\"\n",
    "        super().__init__(A, b)\n",
    "        self.w = w\n",
    "        \n",
    "    def _B_f(self):\n",
    "        D, L, U = self._dlu(self.A)\n",
    "        w = self.w\n",
    "    \n",
    "        _inv_DsWL = np.linalg.pinv(D - w * L)\n",
    "\n",
    "        B = _inv_DsWL @ ((1-w) * D + w * U)\n",
    "        f = w * _inv_DsWL @ self.b\n",
    "    \n",
    "        return B, f\n",
    "    \n",
    "    def method_name(self):\n",
    "        return super().method_name() + f' (w={self.w})'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "JacobiIteration get result x = [-1. -1. -1. -1.] after 68 iterations.\n",
      "GaussSeidel get result x = [-1. -1. -1. -1.] after 36 iterations.\n",
      "SOR (w=0.9) get result x = [-1. -1. -1. -1.] after 46 iterations.\n",
      "SOR (w=1) get result x = [-1. -1. -1. -1.] after 36 iterations.\n",
      "SOR (w=1.1) get result x = [-1. -1. -1. -1.] after 28 iterations.\n"
     ]
    }
   ],
   "source": [
    "A = [[-4, 1, 1, 1], \n",
    "     [1, -4, 1, 1],\n",
    "     [1, 1, -4, 1],\n",
    "     [1, 1, 1, -4]]\n",
    "b = [1, 1, 1, 1]\n",
    "\n",
    "x0 = [0, 0, 0, 0]\n",
    "\n",
    "for method in [JacobiIteration(A, b), GaussSeidel(A, b), SOR(A, b, 0.9), SOR(A, b, 1), SOR(A, b, 1.1)]:\n",
    "    method(x0=x0, eps=1e-9, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "寻找最好的 w："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SOR (w=0.1) get result x = [-0.99999996 -0.99999996 -0.99999996 -0.99999996] after 649 iterations.\n",
      "SOR (w=0.1181818181818182) get result x = [-0.99999997 -0.99999997 -0.99999997 -0.99999997] after 550 iterations.\n",
      "SOR (w=0.13636363636363635) get result x = [-0.99999997 -0.99999997 -0.99999997 -0.99999997] after 476 iterations.\n",
      "SOR (w=0.15454545454545454) get result x = [-0.99999998 -0.99999998 -0.99999998 -0.99999998] after 419 iterations.\n",
      "SOR (w=0.17272727272727273) get result x = [-0.99999998 -0.99999998 -0.99999998 -0.99999998] after 374 iterations.\n",
      "SOR (w=0.19090909090909092) get result x = [-0.99999998 -0.99999998 -0.99999998 -0.99999998] after 337 iterations.\n",
      "SOR (w=0.2090909090909091) get result x = [-0.99999998 -0.99999998 -0.99999998 -0.99999998] after 306 iterations.\n",
      "SOR (w=0.22727272727272727) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 280 iterations.\n",
      "SOR (w=0.24545454545454545) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 258 iterations.\n",
      "SOR (w=0.26363636363636367) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 239 iterations.\n",
      "SOR (w=0.28181818181818186) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 222 iterations.\n",
      "SOR (w=0.3) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 207 iterations.\n",
      "SOR (w=0.3181818181818182) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 194 iterations.\n",
      "SOR (w=0.33636363636363636) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 182 iterations.\n",
      "SOR (w=0.3545454545454545) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 171 iterations.\n",
      "SOR (w=0.3727272727272727) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 162 iterations.\n",
      "SOR (w=0.3909090909090909) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 153 iterations.\n",
      "SOR (w=0.40909090909090906) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 145 iterations.\n",
      "SOR (w=0.42727272727272725) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 137 iterations.\n",
      "SOR (w=0.44545454545454544) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 130 iterations.\n",
      "SOR (w=0.4636363636363636) get result x = [-0.99999999 -0.99999999 -0.99999999 -0.99999999] after 124 iterations.\n",
      "SOR (w=0.4818181818181818) get result x = [-0.99999999 -0.99999999 -0.99999999 -1.        ] after 118 iterations.\n",
      "SOR (w=0.5) get result x = [-1. -1. -1. -1.] after 113 iterations.\n",
      "SOR (w=0.5181818181818182) get result x = [-1. -1. -1. -1.] after 108 iterations.\n",
      "SOR (w=0.5363636363636364) get result x = [-1. -1. -1. -1.] after 103 iterations.\n",
      "SOR (w=0.5545454545454546) get result x = [-1. -1. -1. -1.] after 98 iterations.\n",
      "SOR (w=0.5727272727272728) get result x = [-1. -1. -1. -1.] after 94 iterations.\n",
      "SOR (w=0.5909090909090909) get result x = [-1. -1. -1. -1.] after 90 iterations.\n",
      "SOR (w=0.609090909090909) get result x = [-1. -1. -1. -1.] after 87 iterations.\n",
      "SOR (w=0.6272727272727272) get result x = [-1. -1. -1. -1.] after 83 iterations.\n",
      "SOR (w=0.6454545454545454) get result x = [-1. -1. -1. -1.] after 80 iterations.\n",
      "SOR (w=0.6636363636363636) get result x = [-1. -1. -1. -1.] after 76 iterations.\n",
      "SOR (w=0.6818181818181818) get result x = [-1. -1. -1. -1.] after 73 iterations.\n",
      "SOR (w=0.7) get result x = [-1. -1. -1. -1.] after 71 iterations.\n",
      "SOR (w=0.7181818181818181) get result x = [-1. -1. -1. -1.] after 68 iterations.\n",
      "SOR (w=0.7363636363636363) get result x = [-1. -1. -1. -1.] after 65 iterations.\n",
      "SOR (w=0.7545454545454545) get result x = [-1. -1. -1. -1.] after 63 iterations.\n",
      "SOR (w=0.7727272727272727) get result x = [-1. -1. -1. -1.] after 60 iterations.\n",
      "SOR (w=0.7909090909090909) get result x = [-1. -1. -1. -1.] after 58 iterations.\n",
      "SOR (w=0.8090909090909091) get result x = [-1. -1. -1. -1.] after 56 iterations.\n",
      "SOR (w=0.8272727272727273) get result x = [-1. -1. -1. -1.] after 54 iterations.\n",
      "SOR (w=0.8454545454545453) get result x = [-1. -1. -1. -1.] after 52 iterations.\n",
      "SOR (w=0.8636363636363635) get result x = [-1. -1. -1. -1.] after 50 iterations.\n",
      "SOR (w=0.8818181818181817) get result x = [-1. -1. -1. -1.] after 48 iterations.\n",
      "SOR (w=0.8999999999999999) get result x = [-1. -1. -1. -1.] after 46 iterations.\n",
      "SOR (w=0.9181818181818181) get result x = [-1. -1. -1. -1.] after 44 iterations.\n",
      "SOR (w=0.9363636363636363) get result x = [-1. -1. -1. -1.] after 42 iterations.\n",
      "SOR (w=0.9545454545454545) get result x = [-1. -1. -1. -1.] after 41 iterations.\n",
      "SOR (w=0.9727272727272727) get result x = [-1. -1. -1. -1.] after 39 iterations.\n",
      "SOR (w=0.9909090909090909) get result x = [-1. -1. -1. -1.] after 37 iterations.\n",
      "SOR (w=1.009090909090909) get result x = [-1. -1. -1. -1.] after 36 iterations.\n",
      "SOR (w=1.0272727272727273) get result x = [-1. -1. -1. -1.] after 34 iterations.\n",
      "SOR (w=1.0454545454545454) get result x = [-1. -1. -1. -1.] after 33 iterations.\n",
      "SOR (w=1.0636363636363637) get result x = [-1. -1. -1. -1.] after 31 iterations.\n",
      "SOR (w=1.0818181818181818) get result x = [-1. -1. -1. -1.] after 30 iterations.\n",
      "SOR (w=1.1) get result x = [-1. -1. -1. -1.] after 28 iterations.\n",
      "SOR (w=1.1181818181818182) get result x = [-1. -1. -1. -1.] after 27 iterations.\n",
      "SOR (w=1.1363636363636365) get result x = [-1. -1. -1. -1.] after 25 iterations.\n",
      "SOR (w=1.1545454545454545) get result x = [-1. -1. -1. -1.] after 24 iterations.\n",
      "SOR (w=1.1727272727272728) get result x = [-1. -1. -1. -1.] after 22 iterations.\n",
      "SOR (w=1.190909090909091) get result x = [-1. -1. -1. -1.] after 21 iterations.\n",
      "SOR (w=1.2090909090909092) get result x = [-1. -1. -1. -1.] after 19 iterations.\n",
      "SOR (w=1.2272727272727273) get result x = [-1. -1. -1. -1.] after 17 iterations.\n",
      "SOR (w=1.2454545454545456) get result x = [-1. -1. -1. -1.] after 17 iterations.\n",
      "SOR (w=1.2636363636363637) get result x = [-1. -1. -1. -1.] after 18 iterations.\n",
      "SOR (w=1.2818181818181817) get result x = [-1. -1. -1. -1.] after 19 iterations.\n",
      "SOR (w=1.3) get result x = [-1. -1. -1. -1.] after 20 iterations.\n",
      "SOR (w=1.3181818181818181) get result x = [-1. -1. -1. -1.] after 21 iterations.\n",
      "SOR (w=1.3363636363636364) get result x = [-1. -1. -1. -1.] after 22 iterations.\n",
      "SOR (w=1.3545454545454545) get result x = [-1. -1. -1. -1.] after 23 iterations.\n",
      "SOR (w=1.3727272727272728) get result x = [-1. -1. -1. -1.] after 24 iterations.\n",
      "SOR (w=1.3909090909090909) get result x = [-1. -1. -1. -1.] after 25 iterations.\n",
      "SOR (w=1.4090909090909092) get result x = [-1. -1. -1. -1.] after 26 iterations.\n",
      "SOR (w=1.4272727272727272) get result x = [-1. -1. -1. -1.] after 28 iterations.\n",
      "SOR (w=1.4454545454545455) get result x = [-1. -1. -1. -1.] after 29 iterations.\n",
      "SOR (w=1.4636363636363636) get result x = [-1. -1. -1. -1.] after 31 iterations.\n",
      "SOR (w=1.481818181818182) get result x = [-1. -1. -1. -1.] after 32 iterations.\n",
      "SOR (w=1.5) get result x = [-1. -1. -1. -1.] after 34 iterations.\n",
      "SOR (w=1.5181818181818183) get result x = [-1. -1. -1. -1.] after 35 iterations.\n",
      "SOR (w=1.5363636363636364) get result x = [-1. -1. -1. -1.] after 37 iterations.\n",
      "SOR (w=1.5545454545454547) get result x = [-1. -1. -1. -1.] after 39 iterations.\n",
      "SOR (w=1.5727272727272728) get result x = [-1. -1. -1. -1.] after 41 iterations.\n",
      "SOR (w=1.5909090909090908) get result x = [-1. -1. -1. -1.] after 44 iterations.\n",
      "SOR (w=1.6090909090909091) get result x = [-1. -1. -1. -1.] after 46 iterations.\n",
      "SOR (w=1.6272727272727272) get result x = [-1. -1. -1. -1.] after 50 iterations.\n",
      "SOR (w=1.6454545454545455) get result x = [-1. -1. -1. -1.] after 53 iterations.\n",
      "SOR (w=1.6636363636363636) get result x = [-1. -1. -1. -1.] after 55 iterations.\n",
      "SOR (w=1.6818181818181819) get result x = [-1. -1. -1. -1.] after 60 iterations.\n",
      "SOR (w=1.7) get result x = [-1. -1. -1. -1.] after 64 iterations.\n",
      "SOR (w=1.7181818181818183) get result x = [-1. -1. -1. -1.] after 69 iterations.\n",
      "SOR (w=1.7363636363636363) get result x = [-1. -1. -1. -1.] after 74 iterations.\n",
      "SOR (w=1.7545454545454546) get result x = [-1. -1. -1. -1.] after 81 iterations.\n",
      "SOR (w=1.7727272727272727) get result x = [-1. -1. -1. -1.] after 89 iterations.\n",
      "SOR (w=1.790909090909091) get result x = [-1. -1. -1. -1.] after 97 iterations.\n",
      "SOR (w=1.809090909090909) get result x = [-1. -1. -1. -1.] after 107 iterations.\n",
      "SOR (w=1.8272727272727274) get result x = [-1. -1. -1. -1.] after 120 iterations.\n",
      "SOR (w=1.8454545454545455) get result x = [-1. -1. -1. -1.] after 136 iterations.\n",
      "SOR (w=1.8636363636363635) get result x = [-1. -1. -1. -1.] after 155 iterations.\n",
      "SOR (w=1.8818181818181818) get result x = [-1. -1. -1. -1.] after 181 iterations.\n",
      "SOR (w=1.9) get result x = [-1. -1. -1. -1.] after 216 iterations.\n"
     ]
    }
   ],
   "source": [
    "A = [[-4, 1, 1, 1], \n",
    "     [1, -4, 1, 1],\n",
    "     [1, 1, -4, 1],\n",
    "     [1, 1, 1, -4]]\n",
    "b = [1, 1, 1, 1]\n",
    "\n",
    "for w in np.linspace(0.1, 1.9, 100):\n",
    "    rs = SOR(A, b, w)(x0=[0, 0, 0, 0], eps=1e-9, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "把上面输出的结果写到一个文件 'sorres.txt' 中。然后从中读取 w 与迭代次数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "pattern = r'SOR \\(w=(.*?)\\) get result .*? after (\\d*?) iterations.'\n",
    "\n",
    "ws = []\n",
    "iters = []\n",
    "\n",
    "with open('sorres.txt')  as f:\n",
    "    for line in f:\n",
    "        w, i = list(map(lambda x: float(x), \n",
    "                            re.findall(pattern, line)[0]))\n",
    "        ws.append(w)\n",
    "        iters.append(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作图分析："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHl9JREFUeJzt3X20XXV95/H3B8JTICEkuYQQ8gCaJSJqhVsbK7ZF0UFsDXbUYoMEixOrjtXF6IiTNZ22U9bY6RodXTq0QZwJcAsqlhKVomlgREdAbxSRR4mRQCKQSyIBGoMGvvPH/h2zc3Luvfvee/bZ+5zzea111tn7t/c555udc+/3/h63IgIzM7NmB1UdgJmZ1ZMThJmZteQEYWZmLTlBmJlZS04QZmbWkhOEmZm15ARhZmYtOUGYmVlLThBmZtbStKoDmIq5c+fGkiVLqg7DzKyrbNy48YmIGBjvvK5OEEuWLGF4eLjqMMzMuoqkLUXOcxOTmZm15ARhZmYtOUGYmVlLThBmZtaSE4SZmbXUdwliaAiWLIGDDsqeh4aqjsjMrJ66epjrRA0NwapVsHt3tr9lS7YPsGJFdXGZmdVRX9UgVq/elxwadu/Oys3MbH99lSAefnhi5WZm/ayvEsSiRRMrNzPrZ32VIC69FKZP379s+vSs3MzM9tdXCWLFClizBo4/PtufMyfbdwe1mdmB+ipBQJYMfvzjbPvDH3ZyMDMbTd8lCIAjj4SZM+HRR6uOxMysvvoyQQDMnw8/+1nVUZiZ1VffJojjj3cNwsxsLH2bIFyDMDMbW18niEcfhYiqIzEzq6e+TRDHHw979sCuXVVHYmZWT32bIObPz57dzGRm1lqpCULSLEnXSbpf0n2SXiVptqT1kh5Mz8ekcyXp05I2SbpL0mllxtaYLOeOajOz1squQXwKuCkiTgZeDtwHXAJsiIilwIa0D/BGYGl6rAIuKzMw1yDMzMZWWoKQdDTwO8AVABHxy4h4ElgOrE2nrQXOTdvLgSsjczswS9L8suJrJAjXIMzMWiuzBnEiMAL8b0k/kPQ5SUcC8yKi8Wv5MWBe2l4APJJ7/dZUVooZM+Coo1yDMDMbTZkJYhpwGnBZRLwC+Ff2NScBEBEBTGigqaRVkoYlDY+MjEwpwMZQVzMzO1CZCWIrsDUi7kj715EljMcbTUfpeXs6vg1YmHv9CalsPxGxJiIGI2JwYGBgSgF6NrWZ2ehKSxAR8RjwiKQXpaLXAfcC64CVqWwlcEPaXgdckEYzLQN25ZqiSuHZ1GZmo5tW8vt/ABiSdCiwGXgXWVL6oqSLgC3A29O5NwLnAJuA3encUjVqEBEglf1pZmbdpdQEERF3AoMtDr2uxbkBvL/MeJrNnw+7d8NTT8HRR3fyk83M6q9vZ1KDh7qamY2lrxNEYza1+yHMzA7U1wnCNQgzs9H1dYLwekxmZqPr6wQxYwZMn+4mJjOzVvo6QUieTW1mNpq+ThCQNTO5BmFmdqC+ThBDQ7BxI9x6KyxZku2bmVmmbxPE0BCsWpVNlAPYsiXbd5IwM8v0bYJYvXpfcmjYvTsrNzOzPk4QDz88sXIzs37Ttwli0aKJlZuZ9Zu+TRCXXprNgcibPj0rNzOzPk4QK1bAmjX7agwzZmT7K1ZUG5eZWV30bYKALBls2QIvfSmceaaTg5lZXl8niIaTToLNm6uOwsysXpwg2JcgIqqOxMysPpwgyBLE7t2wfXvVkZiZ1YcTBFmCADczmZnlOUHgBGFm1ooTBNlCfeAEYWaW5wQBHH44LFjgBGFmlldqgpD0kKQfSbpT0nAqmy1pvaQH0/MxqVySPi1pk6S7JJ1WZmzNPNTVzGx/nahBnBkRvxERg2n/EmBDRCwFNqR9gDcCS9NjFXBZB2L7NScIM7P9VdHEtBxYm7bXAufmyq+MzO3ALEnzOxXUSSfBtm2wZ0+nPtHMrN7KThABfEPSRkmrUtm8iGjcBfoxYF7aXgA8knvt1lTWESedlE2U27KlU59oZlZv00p+/zMiYpukY4H1ku7PH4yIkDSh+csp0awCWNTGtbnzQ11f9KK2va2ZWdcqtQYREdvS83bgeuCVwOONpqP03Ji/vA1YmHv5Cams+T3XRMRgRAwODAy0LVbPhTAz219pCULSkZJmNLaBNwB3A+uAlem0lcANaXsdcEEazbQM2JVriirdvHlwxBFOEGZmDWU2Mc0DrpfU+Jx/iIibJH0P+KKki4AtwNvT+TcC5wCbgN3Au0qM7QCSRzKZmeWVliAiYjPw8hblO4DXtSgP4P1lxTOeoaEsOdxzTzaz+tJLfX8IM+tvnklNlhxWrYJf/CLb37Il2x8aqjYuM7MqOUEAq1dny33n7d6dlZuZ9SsnCODhhydWbmbWD5wggNGmU7RxmoWZWddxgiDrkJ4+ff+y6dOzcjOzfuUEQTZaac0aWLw42z/ssGzfo5jMrJ85QSQrVsBDD8G73w0zZsAf/3HVEZmZVcsJosmpp8ITT8D27eOfa2bWy5wgmpx6avZ8993VxmFmVjUniCZOEGZmGSeIJsceC3PnOkGYmTlBNJHgpS91gjAzc4Jo4dRTswTx/PNVR2JmVh0niBZOPRWeecZLbZhZf3OCaMEd1WZmThAtveQl2bMThJn1MyeIFo4+GhYudIIws/7mBDGK2bPhi1+Egw7K7jDnmweZWb8p857UXWtoKLv16N692X7jDnPgBfzMrH+4BtHC6tX7kkOD7zBnZv3GCaIF32HOzMwJoiXfYc7MrAMJQtLBkn4g6atp/0RJd0jaJOkLkg5N5Yel/U3p+JKyYxuN7zBnZtaZGsQHgfty+38DfDIiXgj8HLgolV8E/DyVfzKdV4nGHeZmzsz2Fy70HebMrP+UmiAknQC8Cfhc2hfwWuC6dMpa4Ny0vTztk46/Lp1fiRUr4Mors+1rr3VyMLP+U3YN4n8C/xFoLHs3B3gyIhpjhLYCC9L2AuARgHR8Vzq/Mr/1W9nz7bdXGYWZWTUKJQhJ/13STEmHSNogaUTS+eO85veB7RGxsS2R7nvfVZKGJQ2PjIy0860PcNxx2SQ5Jwgz60dFaxBviIingN8HHgJeCHxknNe8GnizpIeAa8malj4FzJLUmKB3ArAtbW8DFgKk40cDO5rfNCLWRMRgRAwODAwUDH/yli1zgjCz/lQ0QTR+ob8J+FJE7BrvBRHxsYg4ISKWAOcBN0fECuAW4K3ptJXADWl7XdonHb85IqJgfKVZtgweeQS2bRv/XDOzXlI0QXxV0v3A6cAGSQPAnkl+5keBiyVtIutjuCKVXwHMSeUXA5dM8v3batmy7PmOO6qNw8ys01T0j3RJs4FdEfGcpOnAzIh4rNToxjE4OBjDw8Olfsazz2bDXT/0Ifibygbempm1j6SNETE43nkTGcV0MvBHki4gawJ6w2SD6yaHHZbNg/jUp7yyq5n1l0KruUq6CngBcCfwXCoO4MqS4qqNoaFsNVev7Gpm/aZQE5Ok+4BT6tBpnNeJJqYlS7Kk0GzxYnjooVI/2sysFO1uYrobOG5qIXUnr+xqZv2q6A2D5gL3Svou8GyjMCLeXEpUNbJoUesahFd2NbNeVzRB/EWZQdTZpZdmfQ67d+8r88quZtYPCjUxRcQ3gfuBGelxXyrreY2VXefOzfbnzfPKrmbWH4quxfR24LvA24C3A3dIeuvYr+odK1bAAw+ABO97n5ODmfWHop3Uq4HfjIiVEXEB8ErgP5cXVv3Mng2nnw7/8i9VR2Jm1hlFE8RBEbE9t79jAq/tGWedBbfdBs88U3UkZmblK/pL/iZJX5d0oaQLga8BN5YXVj2ddVY2Ye7WW6uOxMysfEU7qT8CrAFelh5rIuKjZQZWR69+NRx+uJuZzKw/FG4miogvR8TF6XF9mUHV1eGHwwteAJ/5jNdlMrPeN+Y8CEnfjogzJD1NtvbSrw8BEREzS42uZoaG4Mc/hl/9Ktv3ukxm1ssKL/ddR51YiynP6zKZWS9o61pMaTXXcct6nddlMrN+UrQP4iX5nXTP6NPbH069jbb+ktdlMrNeNGaCkPSx1P/wMklPpcfTwOPsu5d037j00mwdpjyvy2RmvWrMBBER/y0iZgB/GxEz02NGRMyJiI91KMbaaKzL1KgxHHmk12Uys941kXtSHwMsBQ5vlEVEpVPGOt1JnbdiBXz96/DYYzCt6Jq4ZmY10O5O6ncDtwJfB/4yPf/FVALsdm95C+zYAd/+dtWRmJmVo2gn9QeB3wS2RMSZwCuAJ0uLqgucfXY2ce76vpwyaGb9oGiC2BMRewAkHRYR9wMvKi+s+jvqKDjlFPjsZz2r2sx6U9EEsVXSLOCfgPWSbgBaTBnbR9Lhkr4r6YeS7pH0l6n8REl3SNok6QuSDk3lh6X9Ten4ksn/s8o3NAR33w3PPQcR+2ZVO0mYWa+Y8ExqSb8LHA3cFBG/HOM8AUdGxDOSDgG+TdZUdTHwjxFxraS/A34YEZdJeh/wsoj4U0nnAW+JiD8aK5YqO6k9q9rMulXbOqklHSzp/sZ+RHwzItaNlRzSeRERjTsnHJIeAbwWuC6VrwXOTdvL0z7p+OtSkqklz6o2s143boKIiOeAByRNeL5wSi53AtuB9cBPgCcjYm86ZSuwIG0vAB5Jn7kX2AXMafGeqyQNSxoeGRmZaEht41nVZtbrivZBHAPcI2mDpHWNx3gviojnIuI3gBPIblN68hRibbznmogYjIjBgYGBqb7dpHlWtZn1uqJTvKZ0/+mIeFLSLcCrgFmSpqVawgnAtnTaNmAhWYf4NLJ+jh1T+dwyNWZPr169ry/iE5/wrGoz6x1F7yj3TeAh4JC0/T3g+2O9RtJAGvmEpCOA1wP3AbcAb02nrWTfmk7r0j7p+M1R87XIV6zIOqTvvDPbb9wnwsysFxSdSf3vyDqO/z4VLSAb8jqW+cAtku4iSyjrI+KrwEeBiyVtIutjuCKdfwUwJ5VfDFwykX9IlV7+8mz00sUXe06EmfWOok1M7yfrQ7gDICIelHTsWC+IiLvIZlw3l29O79Vcvgd4W8F4amVoCH72M99pzsx6S9FO6mfzw1pTH0Gtm386afXqA5uXdu/Oys3MulXRBPFNSf8JOELS64EvAV8pL6zu4jkRZtaLiiaIS4AR4EfAe4AbI8J/HyeeE2FmvahogvhARFweEW+LiLdGxOWSPlhqZF3EcyLMrBcVTRArW5Rd2MY4ulrjTnOLF0NjcZALL3QHtZl1t/HuSf0OSV8BTszPoE6T3nZ2JsTu0JgTsXcvHHccXH65h7yaWXcbb5jrd4BHgbnA/8iVPw3cVVZQ3eyaa2DnTg95NbPuN+HlvuukyuW+R+NlwM2s7oou9z1mDULStyPiDElPs/+8B5Gt6D1zinH2HA95NbNeMWaCiIgz0vOMzoTT/RYtal2D8JBXM+s2RUcxWUGthrxCljTcYW1m3cQJos3yQ16b+b7VZtZNnCBK0Bjy2ipJeI0mM+sWThAlcoe1mXUzJ4gSeY0mM+tmThAlatVhfcQRXqPJzLqDE0SJWq3RFAHvfKdHNJlZ/TlBlKzRYX3VVTBtGuzZkyUJj2gys7pzguiQ1auzhfzyPKLJzCZqaChrgejEYqBF70ltU+QRTWY2VUNDWcvD7t3ZftmLgboG0SGjjVyKcH+EmRWzevW+5NBQZkuEE0SHjLYEB7g/wsyK6XRLRGkJQtJCSbdIulfSPY1blEqaLWm9pAfT8zGpXJI+LWmTpLsknVZWbFUYawkOcH+EmY2v03OryqxB7AX+Q0ScAiwD3i/pFOASYENELAU2pH2ANwJL02MVcFmJsVWiMaKpMeS1mfsjzGwsf/3XB5ZNn17e3KrSEkREPBoR30/bTwP3AQuA5cDadNpa4Ny0vRy4MjK3A7MkzS8rvip5hrWZTcaLX5w9z5mT/aG5eHHWMlHW3So70gchaQnwCuAOYF5EPJoOPQbMS9sLgEdyL9uaynqOlwQ3s8m46abs+d574fnnsxaJMm9lXHqCkHQU8GXgQxHxVP5YZPc7ndA9TyWtkjQsaXhkZKSNkXaOlwQ3s8m46SY47TQ49tjOfF6pCULSIWTJYSgi/jEVP95oOkrP21P5NmBh7uUnpLL9RMSaiBiMiMGBgYHygi+ZlwQ3s4nYtQtuuw3OPrtzn1nmKCYBVwD3RcQncofWASvT9krghlz5BWk00zJgV64pqmeN1jHt5iYzy9uwAZ57rkcSBPBq4J3AayXdmR7nAB8HXi/pQeCstA9wI7AZ2ARcDryvxNhqY6yOaTc3mRlkvwNWpj+rzz+/c78TlHUDdKfBwcEYHh6uOowpaZ4638rixVlzlJn1n1a/I6ZPn9roJUkbI2JwvPM8k7pi402gA8+PMOtnnV5eI88JogbG6rAGr9dk1s+qXOjTCaJGvF6TmTVbuLB1eScm1jpB1IjXazKzZn/4hweWlbm8Rp4TRM2Mt16Th7+a9Zef/hRmzsxqDJ1YXiPPCaKmPPzVzLZvh699Dd7znuznvhPLa+Q5QdTUWP0RkDU3nX++axNmvWpoCE4+ObtV8dVXV/Nz7gRRU0WGv4JrE2a9qDH34ec/z/YffbSan3NPlOsCS5ZkiWAsnkxn1jtG+5lv18+5J8r1kPGamyD7Mh10kJuczHpBlXMf8pwgukDR5qYINzmZ9YLRlvPu9E3FnCC6RGP469VXj1+bcAe2WfeKgFmzDizv1NyHPCeILpOvTYw2V6LBtQmz7jE0lP1Rd/DB8MAD8JrX7Ps57+Tchzx3Unc5d2Cbdb8yVmwdizup+0TRDmw3N5nVV5Urto7FCaLLTWS+xDvfmVVXnSzM6qUuo5aaOUH0gKId2I3WRPdNmNXLaKOTOj1qqZkTRA8pWpsAj3Qyq4NGx3SrfsQqRi01c4LoMePdfKiZaxNm1Wh0TOeTQ2NkYlWjlpo5QfSoIp3XDa5NmHVeq47piH2jDqtODuAE0bOam5vGmzMB7sg266S6dkznOUH0sEZzUwRcdVWxZid3ZJuVq9HvMNoUtKo7pvNKSxCSPi9pu6S7c2WzJa2X9GB6PiaVS9KnJW2SdJek08qKq19NZKmOBjc9mbVXq36HvDp0TOeVWYP4P8DZTWWXABsiYimwIe0DvBFYmh6rgMtKjKuvTWSkU4NrE2bt0arfoaEuHdN5pSWIiLgV2NlUvBxYm7bXAufmyq+MzO3ALEnzy4qt37k2YdZZYw1nhazfry4d03md7oOYFxGPpu3HgHlpewHwSO68ranMSuSObLPyjdesBPXqd8irrJM6slUCJ7xSoKRVkoYlDY+MjJQQWX+Zake2k4VZa41aw/nnj96sBPXrd8jrdIJ4vNF0lJ63p/JtwMLceSeksgNExJqIGIyIwYGBgVKD7TeTaXpysjA7UJFaA9Sz3yGv0wliHbAyba8EbsiVX5BGMy0DduWaoqzDJtORDU4WZg1jdUY31GlC3GjKHOZ6DXAb8CJJWyVdBHwceL2kB4Gz0j7AjcBmYBNwOfC+suKyYiZTm8jLJ4t3vQvmzvU9s633jdcZ3VDnZqU83zDIxjU0lP1FtGVLViuY6lem8R6LF2c/JHX+C8qsqFY3/WmlDt973zDI2ma0juwio55acVOU9ZKJdEZffXX9m5XynCBsQpwszPbplc7o0ThB2KQ5WVi/KlprgO7ojB6NE4S1RatkIcGcOXDooRN/PycLq5tGUpCy7+R4tQbons7o0ThBWNs1ksXzz8MTT8DnPz+12kWrZDF3rkdGWflGSwpFBmp0a7NSnhOEla6dTVGNH8wdO7JHhGsZVo7m/oWio/e6sTN6NE4Q1lHt7rdocC3D2mUi/QvNeqHWkOcEYZUpO1m4lmFFTaZ/Ia+Xag15ThBWC2UlizzXMixvKv0LsO+72Wu1hjwnCKudTiaL5lpGflkQJ4/u10gCzf+fc+fCn/zJ1JLCVVdlr+u1WkOeE4TV2ljDZ+fMyc5pZ+L41a/2JYzRmqicOOqtVc2g+f9zxw745S8n9r79khTynCCsazQPn33iiXJrGc3ct1FfU20uGkuv9i8U4QRhXa/TtYxm7tvonFZNRmUkhX7oXyjCCcJ6St1rGU4cxYyXCPJNRtD+pNBvTUmjcYKwvjBeLWMqy4IUMZHE0etJZKyO404kgoZDDtn3f++k0JrvB2GW07j3xcMPw+zZWdmOHe25D8Zk5e+fcc45cOONWXyLFlV/X4Gx1OVaNj6v0dy4c2f9r13Zit4PwgnCrIC6/LJrVsUvv1bXYufOA7effnriI4XaxTelGptvGGTWRlX3bYym3X0ek23+accw0qlyH0L7uQZh1iZ1rWXktapx1DXW0bjJaOpcgzDrsPFqGZ0cejuaVjWOsjqBp6pxffIDCfK1g8Y1fv551xTK4gRhVrKJJI46JJFOaPXL34mgfpwgzCrSKnGMlUQWL4b3vrfaPo+JyA8jLfLL34mgfqZVHYCZjW7FitF/UVbZ5zFWX4b7BHpHrWoQks6W9ICkTZIuqToeszorq89jss0/rgH0ntokCEkHA58F3gicArxD0inVRmXWfSbb5+Ff/tasTk1MrwQ2RcRmAEnXAsuBeyuNyqxHjNVcZdZKbWoQwALgkdz+1lS2H0mrJA1LGh4ZGelYcGZm/aZOCaKQiFgTEYMRMTgwMFB1OGZmPatOCWIbsDC3f0IqMzOzCtQpQXwPWCrpREmHAucB6yqOycysb9Wmkzoi9kr698DXgYOBz0fEPRWHZWbWt7p6sT5JI8CWquMYx1zgiaqDKMBxtle3xAndE6vjbJ/FETFuJ25XJ4huIGm4yKqJVXOc7dUtcUL3xOo4O69OfRBmZlYjThBmZtaSE0T51lQdQEGOs726JU7onlgdZ4e5D8LMzFpyDcLMzFpygpik8ZYml3SxpHsl3SVpg6TFuWPPSbozPUqdDFggzgsljeTieXfu2EpJD6bHyjLjLBjrJ3Nx/ljSk7ljnbymn5e0XdLdoxyXpE+nf8ddkk7LHevYNS0Q54oU348kfUfSy3PHHkrld0oq9cbvBeL8PUm7cv+/f5471rFbBBSI8yO5GO9O38nZ6VjHrmdbRYQfE3yQTeT7CXAScCjwQ+CUpnPOBKan7fcCX8gde6ZGcV4IfKbFa2cDm9PzMWn7mCpjbTr/A2STKTt6TdNn/Q5wGnD3KMfPAf4ZELAMuKOiazpenL/d+HyyZfbvyB17CJhbk+v5e8BXp/qdKTvOpnP/ALi5iuvZzodrEJPz66XJI+KXQGNp8l+LiFsiYnfavZ1sbalOGzfOMfwbYH1E7IyInwPrgbNLihMmHus7gGtKjGdUEXErsHOMU5YDV0bmdmCWpPl0+JqOF2dEfCfFAdV9R4tcz9FM5fs9YROMs7LvZzs5QUxOoaXJcy4i+4uy4fC0ZPntks4tI8CkaJz/NjU1XCepsWDiRP+NU1X481Jz3YnAzbniTl3TIkb7t3T6mk5E83c0gG9I2ihpVUUx5b1K0g8l/bOkl6SyWl5PSdPJEv+Xc8V1u56F1GYtpl4l6XxgEPjdXPHiiNgm6STgZkk/ioifVBMhXwGuiYhnJb0HWAu8tqJYijoPuC4insuV1emadhVJZ5IliDNyxWek63kssF7S/ekv6Cp8n+z/9xlJ5wD/BCytKJYi/gD4fxGRr23U6XoW5hrE5BRamlzSWcBq4M0R8WyjPCK2pefNwP8FXlFVnBGxIxfb54DTi762zSbyeefRVH3v4DUtYrR/S+2WtJf0MrL/9+URsaNRnrue24HryZpzKhERT0XEM2n7RuAQSXOp4fVMxvp+Vn49J6TqTpBufJDVvDaTNXM0Osde0nTOK8g60JY2lR8DHJa25wIPUlLHWsE45+e23wLcnrZnAz9N8R6TtmdXeU3TeSeTdfipimua+8wljN6p+ib276T+bhXXtECci4BNwG83lR8JzMhtfwc4u8I4j2PfnK1XAg+na1voO9OpONPxo8n6KY6s8nq26+EmpkmIUZYml/RXwHBErAP+FjgK+JIkgIcj4s3Ai4G/l/Q8WQ3u4xFRyn23C8b5Z5LeDOwl+2JfmF67U9J/JbtPB8Bfxf5V5ipiheyvs2sj/bQlHbumAJKuIRtZM1fSVuC/AIekf8ffATeSjWTaBOwG3pWOdfSaFojzz4E5wP9K39G9kS0yNw+4PpVNA/4hIm6qMM63Au+VtBf4BXBe+v/v6C0CCsQJ2R9Z34iIf829tKPXs508k9rMzFpyH4SZmbXkBGFmZi05QZiZWUtOEGZm1pIThJmZteQEYWZmLTlBmJlZS04QZm2S7gfwZ2n7k5JuTtuvlTRUbXRmE+cEYdY+3wJek7YHgaMkHZLKar8wm1kzJwiz9tkInC5pJvAscBtZongNWfIw6ypei8msTSLiV5J+Srae1XeAu8juLPhC4L4KQzObFNcgzNrrW8CHyZqUvgX8KfCD8KJn1oWcIMza61vAfOC2iHgc2IObl6xLeTVXMzNryTUIMzNryQnCzMxacoIwM7OWnCDMzKwlJwgzM2vJCcLMzFpygjAzs5acIMzMrKX/D4ALfewQesyYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(ws, iters, 'bo-')\n",
    "plt.xlabel('w')\n",
    "plt.ylabel('iterations')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从中找到使迭代次数最少的 w:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best w: 1.2272727272727273\n",
      "iterations: 17\n"
     ]
    }
   ],
   "source": [
    "w_best,iters_min = min(zip(ws, iters), key=lambda x: x[1])\n",
    "print(f'best w: {w_best}\\niterations: {int(iters_min)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一个无法直接解的例子："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  2., -1.])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[1, -3, -6], \n",
    "     [2, 8, -3], \n",
    "#      [5, 2, 1]]\n",
    "b = [1, 21, 8]\n",
    "np.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:48: RuntimeWarning: overflow encountered in matmul\n",
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:48: RuntimeWarning: invalid value encountered in matmul\n",
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:50: RuntimeWarning: invalid value encountered in less_equal\n"
     ]
    },
    {
     "ename": "Exception",
     "evalue": "cannot reach eps (1e-08) after max_steps (5000). The last result: x = [nan nan nan]",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-53-1c82c63562b5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgauss_seidel_iter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-8\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-45-6fde854fad32>\u001b[0m in \u001b[0;36mgauss_seidel_iter\u001b[0;34m(A, b, x0, eps, max_steps, verbose)\u001b[0m\n\u001b[1;32m     51\u001b[0m             \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     52\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     55\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: cannot reach eps (1e-08) after max_steps (5000). The last result: x = [nan nan nan]"
     ]
    }
   ],
   "source": [
    "gauss_seidel_iter(A, b, eps=1e-8, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "做一些行变化之后就可以解了："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gauss_seidel_iter get result x = [ 1.  2. -1.] after 12 iterations.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 1.,  2., -1.])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A0 = [[5, 2, 1], \n",
    "      [2, 8, -3], \n",
    "      [1, -3, -6]]\n",
    "b0 = [8, 21, 1]\n",
    "gauss_seidel_iter(A1, b1, eps=1e-8, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "比较 J 和 GS 的收敛性："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = [[1, 2, -2], [1, 1, 1], [2, 2, 1]]\n",
    "\n",
    "A2 = [[2, -1, 1], [2, 2, 2], [-1, -1, 2]]\n",
    "\n",
    "b = [1, 1, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3.,  3.,  1.])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.solve(A1, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.16666667, -0.16666667,  0.5       ])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.solve(A2, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jacobi: jacobi_iter get result x = [-3.  3.  1.] after 3 iterations.\n",
      "gauss_seidel: "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:48: RuntimeWarning: overflow encountered in matmul\n",
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:48: RuntimeWarning: invalid value encountered in matmul\n",
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:50: RuntimeWarning: invalid value encountered in less_equal\n"
     ]
    },
    {
     "ename": "Exception",
     "evalue": "cannot reach eps (1e-05) after max_steps (5000). The last result: x = [nan nan nan]",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-65-a8fd089a47a0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjacobi_iter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'gauss_seidel: '\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgauss_seidel_iter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-45-6fde854fad32>\u001b[0m in \u001b[0;36mgauss_seidel_iter\u001b[0;34m(A, b, x0, eps, max_steps, verbose)\u001b[0m\n\u001b[1;32m     51\u001b[0m             \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     52\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     55\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: cannot reach eps (1e-05) after max_steps (5000). The last result: x = [nan nan nan]"
     ]
    }
   ],
   "source": [
    "print('jacobi: ', end='')\n",
    "j = jacobi_iter(A1, b, verbose=True)\n",
    "print('gauss_seidel: ', end='')\n",
    "g = gauss_seidel_iter(A1, b, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gauss_seidel: gauss_seidel_iter get result x = [ 0.16666863 -0.16666896  0.49999983] after 19 iterations.\n",
      "jacobi: "
     ]
    },
    {
     "ename": "Exception",
     "evalue": "cannot reach eps (1e-05) after max_steps (5000). The last result: x = [7.99615552e+241 1.08514261e+242 2.57044248e+241]",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-67-15461764c0a3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgauss_seidel_iter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'jacobi: '\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjacobi_iter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-4-d9fdc2d22bab>\u001b[0m in \u001b[0;36mjacobi_iter\u001b[0;34m(A, b, x0, eps, max_steps, verbose)\u001b[0m\n\u001b[1;32m     52\u001b[0m             \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     53\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 54\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"cannot reach eps ({eps}) after max_steps ({max_steps}). The last result: x = {x}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     55\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     56\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: cannot reach eps (1e-05) after max_steps (5000). The last result: x = [7.99615552e+241 1.08514261e+242 2.57044248e+241]"
     ]
    }
   ],
   "source": [
    "print('gauss_seidel: ', end='')\n",
    "# g = gauss_seidel_iter(A2, b, verbose=True)\n",
    "print('jacobi: ', end='')\n",
    "j = jacobi_iter(A2, b, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
